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Abstract 

We have developed a couple of optimal damping algorithms (ODAs) for unrestricted 
Hartree-Fock (UHF) calculations of open-shell molecular systems. A series of equations 
were derived for both concurrent and alternate constructions of a and /3 Fock matrices 
in the integral-direct self-consistent-field (SCF) procedure. Several test calculations were 
performed to check the convergence behaviors. It was shown that the concurrent algorithm 
provides better performance than does the alternate one. 
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1. Introduction 

The Hartree-Fock (HF) method has long been the fundament in molecular orbital 
(MO) calculations. In 1951, Roothaan [1] developed systematic operator formulations 
known as the restricted HF (RHF) method for the ground state of closed-shell systems 
where each occupied orbital has a pair of a- and /3-spin electrons. To treat radicals with 
open-shells, Pople and Nesbet [2] provided the unrestricted extension (UHF) by which 
CK-spin orbitals and /3-spin orbitals can be different by incorporating the spin polarization 
at the cost of spin contamination. Roothaan [3] successively proposed the restricted ver- 
sion for open-shell molecules (ROHF) without the contamination problem, and afterward 
Plakhutin et al. [4] remade the operator formahsm of ROHF in a canonical style. Re- 
cently, the close relationship between UHF and ROHF was re-evaluated by Tuchimochi 
and Scuseria [5,6]. 

The HF equations should be solved under the self-consistent-field (SCF) condition with 
nonlinear dependence through the density matrix. However, it is well recognized that the 
simple iterations suffer from the convergence difficulty even for closed-shell cases [7-9]. 
To reduce the difficulty, there were two older damping techniques proposed by Saunders 
and Hiller [10] (level shift) and by Zerner and Hehenberger [11] (dynamic damping). 
Then, Pulay [12,13] invented a breakthrough approach to improve the convergence of HF 
SCF procedure, named as the direct inversion in the iterative subspace (DIIS). DIIS was 
designed to minimize the squared norm of residuum under a normalization constraint. 
After Pulay's success, various related variants were developed [14-21], where these DIIS 
methods were usually formulated for the Fock and density matrices with atomic orbital 
(AO) basis functions to expand MOs. Particularly, C2-DIIS by Sellers [15] as well as 
Energy-DIIS (EDIIS) by Kudin et al. [16] have been used most widely. The latter 
has a connection to the optimal damping algorithm (ODA) derived by Cances and Bris 
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[22,23], a direct minimization technique of RHF energy under the relaxed constraints of 
idempotency. 

In the direction of MO-oriented SCF optimization, Bacskay [24] pioneered the quadrat- 
ically convergent SCF (QC-SCF) procedure through the second-order Newton optimiza- 
tion of the RHF energy. In QC-SCF, the occupied MOs should be improved by the explicit 
mixing with the unoccupied MOs without any Fock matrix diagonalization with respect 
to AO-indices during the iteration. However, the actual computations could be too costly 
due to the explicit evaluation of orbital Hessian. Anyhow, the diagonalization-free nature 
of MO optimization is favorable to the integral-direct SCF calculation [25] with paral- 
lelism [26]. Hence, several less expensive procedures with efficient approximations [27-31] 
have thus been devised. An AO-based Newton technique [32] was developed as well. 

As pointed out in Ref. [7], UHF could easily encounter the convergence problem. 
Claxton and Smith [33] reported the direct minimization recipe for improvement. Seeger 
and Pople [34] proposed an MO-based optimization approach for UHF, while Bacskay 
[35] extended his QC-SCF [24] for the ROHF case. Neese [36] revised the approximated 
second-order SCF (SO-SCF) method of Ref. [30] by adding the a-^ coupling elements. 
In comparison with the case of closed-shells or RHF, works oriented to the calculations 
of open-shell (UHF or ROHF) have been rather limited. 

In this paper, we propose a couple of ODAs designed for UHF calculations [2,9]; UHF 
could still be a reasonable zeroth-order treatment for open-shell systems as long as the 
spin contamination is small enough. Cances and Bris [23] certainly addressed the applica- 
tion of their ODA/RHF to the UHF case. Nonetheless, both corresponding formulation 
and numerical result were not shown. We present here detailed formulations and algo- 
rithmic descriptions of ODA/UHF, which should be useful for further methodological 
developments. Our proposal covers two ways of AO-based Fock matrix construction in 
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the integral-direct SCF procedure [25], by incorporating several a-^ coupling terms. One 
is the alternate construction in which Fock matrix for each spin (say a) is computed at 
a certain step of iteration and then the corresponding density matrix is updated for the 
construction of another spin The other is the concurrent construction where both a 
and j3 Fock matrices are simultaneously computed and a and /3 density matrices are then 
updated in a single step; the cost of direct integral generations is a half of that for the 
alternate construction. An attractive point of ODA is a relatively small requirement of 
memory resource. The remaining part of this paper is organized as follows. In Section 2, 
a brief summary of the ODA/RHF method in Ref. [23] is given for self-completeness and 
later convenience. Section 3 describes two ODAs for UHF in detail. In Section 4, test 
applications with four examples, e.g. CN radical, are shown. 

2. Brief summary of ODA/RHF 

Since the basic formulation and notation of our ODA/UHF follow the original ones 
of ODA/RHF by Cances and Bris [23], the essential equations are summarized in this 
SubSection. For simplicity, all matrices are written in Capital italic font (or without 
bold font) hereafter. The subscript specifies the step number of SCF iteration. The total 
number of electrons is 2A^e in the closed-shell RHF description. The number of AOs (or 
dimension of matrix) is N which specifies the formal dimension of matrices. 

The Fock matrix F is given as [9] 

F = h + G{D), (1) 
G{D)^2J{D)-K{D), (2) 
D = CC*. (3) 

Here, h is the one-electron contribution (kinetic energy and nuclear attraction energy), 
and the two-electron contribution G (consisting of Coulomb J and exchange K) has 
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the dependence on the density matrix D formed from the AO-MO coefficient matrix C; 
asterisk corresponds to transposition. C is obtained by solving the general eigenequation 

FC = SCe, (4) 

where S is the overlap matrix and s contains the orbital energies in the diagonal elements. 
D of Eq. (3) satisfies two crucial constraints associated with the orthonormal condition 

Ne^TilDS] (Number of electron), (5) 
DSD^D (Idempotency). (6) 

The RHF electronic energy is then written as 

E^^^{D) = Tr [2hD + G{D)D\ = Tr [/iL> + F{D)D\ . (7) 

C is updated through the iterative SCF optimization until the convergence criteria in- 
volving energy and density are satisfied under the given thresholds. 
The k-ih step of SCF iteration consists of 

1. Assemble Fock matrix Fk — h + G{Dk), 

2. Compute RHF energy Ek^Tr [hD^ + FkDk], 

3. Solve eigenvalue problem F^Ck — SCkSk, 

4. Form density matrix -Dfe+i — CkC^, 

5. Check convergence; go to A; -|- 1-th step if not converged. 

As already denoted in Section 1, the above hsted procedures are generally slow to converge 
[7-9], and thus a variety of acceleration techniques are highly necessary. 
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Now, the ODA/RHF procedure is briefed according to Ref. [23]. Two types of density 
matrix D and D are considered, corresponding to strict and relaxed constraints, respec- 
tively. It is notable that that D satisfies 



A^e = Tr 



DS 



(8) 



like Eq. (5) but that the idempotency requirement is relaxed 



DSD < D 



(9) 



unlike Eq. (6). In ODA, -Dfc+i is defined with an interpolation parameter A G [0, 1] as 



bk+i = (1 - \)Dk + = Dk + \{Dk+i - Dk) = Dk + XADk+i. 



So, the RHF Fock matrix F{D) = h + G{D) of Eq. (1) is deformed as 



(10) 



Fk+i — F{Dk+i) 

= h + (l-X)GiDk) + XG{Dk+i) = {l-X)Fk + XFk+i. (11) 

The parameter A is given by minimizing the RHF energy of Eq. (7). This is simply solved 
by a line search for minimizer 



£;^"'^(i5fe+i) = E^'^' iD,) + sX + cA^ A e [0, 1] 



nRHF 



(12) 



where the parameters s and c are given by 



s = Tr 

c = Tr 



F{Dk)ADk+i 
{F{Dk+i)-F{Dk)}AD 

k+l 



(13) 
(14) 



It is easy to find the optimial value of A by a condition 



dX 



= s + 2cA = 0, 



(15) 



where s should be negative until converge in a sense of steepest descendent slope, as 
carefully discussed in Ref. [29]. Note that c should be positive conversely. The optimal 
damping parameter A e [0, 1] can then be given by 



I —s/2c, otherwise. ^ ^ 

It is notable that ODA does not work in the case of unity in Eq. (16) (see also Eq. (10)). 
As a whole, the ODA-based RHF calculations are sketched as below. 

1. Initialization: Choose an initial guess Dq. Assemble Fq = F{Dq). Compute Eq — 
E{Do). Set Do = Dq, Fq = Fq and A; = 0. 

2. Iteration: 

(a) Diagonalize F^ and assemble the density matrix D^+i via the aufbau principle 



(b) If the difference Dk+i — is enough small then go to termination; the energy 
convergence should be checked as well. 

(c) Assemble the Fock matrix F^+i = F{Dk+i) and compute the RHF energy 



A 





[9]. 



Fk+i^h + G{Dk+i), 



Ek+i — Tr [hDk+i + Fk+iDk+i] . 



(d) 



Set ADk+i = -Dfc - -Dfe+i. 



Compute 



c = Tr [{Ffc+i - Fk} AD, 



k+l 



8 



(f) Set A = lifc<— s/2 and A = —s/2c otherwise, and interpolate 

Dk+i = (1- A)L>fe + AL'fe+i, 
Fk+i = (1 — A)Fjk + AFfe+i, 

(g) Set k = k + 1 and go to 2(a). 

3. Termination: Set C — C^+i, D — Dk+i, F — F^+i and E — Ek+i. 

Once the early stage of iteration with a reasonable guess for density matrix has passed 
successfully, the SCF procedure with ODA [23] can generally be switched to that with one 
of DIIS methods [12,13,15,16], which may be more efficient in accelerations; nonetheless 
ODA is usable for the final convergence, as well. 

3. Proposal of ODA/UHF 

3.1. Alternate version 

The alternate version of UHF/ODA can be regarded as a straightforward extension of 
ODA/RHF [23], except for some points addressed later. The superscripts oi a or (3 for 
matrices specify the spin components, according to Ref . [9] . 
3.1.1. Basic formulation 

First, the numbers of a-spin and /3-spin electrons are set as and A^f , respectively, 
for the density matrices. The a and (3 Fock matrices are defined, respectively, as [9] 

F''^h + G'{D") + J{Df^), (17) 
^h + G'{D^) + J{D''), (18) 
G'{D) = J{D)-K{D), (19) 

where the two density matrices are given in the same way as Eq. (3) 

= C'^C'a*^ (20) 
= C^C^*. (21) 
9 



The AO-MO coefficients are obtained by solving the pair of eigenequations 



The UHF electronic energy is then written as 

^UHFpa^ L)/') = Tr + D^) + ]-G'{D'^)D" + ]-G'{D^)D^ + J{D^)D'^ 



(22) 
(23) 



= iTr 
2 



/i(D" + D^) + F{D")D" + F(D^)D^] , 



(24) 



where an equivalence relation due to a classical nature of Coulomb interaction 



Tr [j(L)^)L)"] = Tr [j{D'^)D 



(25) 



is utilized. 

For a while, we focus on updating the /3 density matrix, assuming that the SCF 
iteration starts on the a Fock matrix construction. When the relaxed /3 density matrix 
incorporating is defined as 



4^^, = (1 - A^)^f + A/'^, = + A/^A<„ 



the resulting (3 Fock matrix with the available strict I^^+i is given by 



(26) 



k+l 



= F(L>f+i,L>?+i), 

= h+{l- A^)G'(^f ) + A^G'«,) + J(D,\,), 
= (l-A^)[/i + G'(Df) + J(D,"); 

+A^ [/i + G'(i)f^J + J{D-,^,) 

+(1 - A^) [J(D,\J - J{Dt) . 
= (l-A/^)Ff + A^<, 

+(1 - A^) [j{Dt^,) - J{Dt) . 
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(27) 
(28) 



(29) 



(30) 



and the ^ Coulomb part associated with the next a Fock matrix becomes 

^f-M = = (1 - A^) Jf + 

The UHF energy minimizer of interest is then defined as 

^^""(^fe+i, Dl.) = ^"""(^fe+i, Di) + /A^ + c^(A^)^ e [0, 1]. 
The crucial parameters and are here given by 



(31) 



(32) 



^Tr [{F(Df+„ Dt^,) - F{bl Dt^,)]AD 



2 
2 



/3 



/3 ■ 
fe+1 



(33) 



(34) 



with the relation 



n^f,^fe\i) = ^fe' + ^fe"+i-^fe"- 
As in the case of ODA/RHF described in the previous SubSection, the condition 



(35) 



2d^>P = 0, 



(36) 



leads to an optimal A'^ which minimizes the UHF energy function of Eq. (32), and the 
result is shown as 



A^ 



1, if \c^\ < \s^\/2 or s^c^ > 

—s^/2c^, otherwise. 



(37) 



The conditions for and are shghtly modified from those of ODA/RHF [29], since the 
simple assumption of a steepest descendent search for the RHF energy is not valid for the 
alternate UHF calculations due to the stepwise coupling with the a matrices. 
3.1.2. Algorithmic flow 

The above-mentioned way to derive A'^ is apphcable to A" as well. The ODA/UHF 
procedure for alternate Fock matrix constructions can now be stated as follows. 
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Initialization: Choose initial guesses and = D^_^. Assemble Fg' = F{D^, D^), 
= F{D^,D^), = J{D^) and = J{D^). Compute Eq = E{D^,D^_^). Set 
= F- = F-, b^^ = dI = AL>o^ = and = 0. 

Iteration: 

(a) Diagonalize F^ and assemble the density matrix D^_^^ via the aufbau principle. 

(b) Assemble the Fock matrix Ff and the Coulomb integrals J^+i- 

Fi = h + G'{D',) + J{Dt^,). 

(c) Set AD^^, = Dt- D'^_,,. 

(d) Compute 

(e) Set A'^ = 1 if |c^| < |s''|/2 or s^c^ > and = -s^ /2c^ otherwise, and 
interpolate 

Fi = (l-AOFf_, + A^Ff + (l-A/^)(J,V-Jn- 

(f ) Diagonalize and assemble the density matrix via aufbau principle. 

(g) Assemble the Fock matrix F^_^-^ and the Coulomb integrals J^+\- 

Fk+i = h + G\Dl^,) + J{Dl^,). 
12 



(h) Set d.Dl^, = Di - Di^,. 

(i) Compute the UHF energy 



(j) If the differences -Dfc+i — -Dfc and -Df+i — -Of are enough small then go to 
termination; the energy convergence should be checked as well. 

(k) Compute 

(1) Set = 1 if |c°| < |s°|/2 or s"c°^ > and A" = -s'^/2c'^ otherwise, and 
interpolate 

DU, = (1 - A")4" + A"D,V, 

^.+1 = (l-A")F," + A"F,V + (l-A°)(jf+,-jf). 

(m) Set k = k + 1 and go to 2(a). 

3. Termination: Set = C^+„ = Cf^^ = D^^„ = D^^,, = F^^,, 
F^ = Ff+i and E = Ek+^. 

As can be seen above, the computational cost of the alternate version of ODA/UHF is 
roughly twice that of ODA/RHF when the integral-direct processing [25] is pursued. 

3.2. Concurrent version 

In the concurrent UHF calculation, both and are simultaneously updated in a 
certain step of SCF iterations. The concurrent ODA/UHF procedure is rather comphcated 
in comparison with the alternate ODA/UHF just shown. This complexity is attributed 
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to the two dimensional nature of minimization algorithm instead of one dimensional line 
search algorithm in the case of alternate version. 

3.2.1. Bcisic formulation 

A couple of Fock matrices including both relaxed a and (3 density matrices are defined 
as (refer also to Eq. (27) for the alternate case) 



~/3 



(38) 
(39) 



F = F{D(^,D^) = h + G'{D^) + JiD""), 
and the updates of density matrices in a step of iteration are done as (see Eq. (26)) 

Dt+i = (1 - A")£'," + A°L>«+, = 4" + A«AL>^+„ A« e [0, 1], (40) 
= (1 - >^')D'k + A^^f+i = + A^A<„ e [0, 1]. (41) 
The a Fock matrix is then calculated as 



F 



k+l 



= F{Dt^„D^,^,), (42) 

= h + {l-X'^)G'{Dt) + X'^G'{Dt^,) + {l-X^)J{D^,) + X^J{D^,^,), (43) 
= (1 - A") [h + G\bt) + J(^f )] + A" [/i + G\Dt^,) + J(Df+,)] 

+(A'^ - A") [j(Df+,) - J(Df )] , (44) 



(1 - A")F, + A'^F.V + (A^ - A") 



f - ~f 



(45) 



and the final expression of the fi Fock matrix becomes 

= (1 - A^)Ff + A^<, + (A-^ - A^) 
The Coulomb matrices are derived as (see Eq. (31)) 



ja ja 



(46) 
(47) 



ja 



f/3 



J(4V) = (1-A")J^ + A'^J,V, 
^(^f+i) = (l-A'^)jf + A'^jf+i. 



(48) 
(49) 
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The UHF energy minimizer is set with an a-^ couphng term as 
where the five crucial parameters s", s^, c", c'^ and i are obtained as 



Tr 
Tr 

2 

^Tr [G(A<,)A<, 



(51) 
(52) 



2 








1 

2 ' 


(53) 


2 








1 

-2^' 


(54) 



i = lY [j(ADf^.,)AD^+, 
= Tr [(Jf-J,^ 



Tr 



AZ^f^, 



(55) 



In the above derivation, an equivalence relation of Eq. (25) is used. The coupling via t 
should be effective in accelerating the SCF convergence, as pointed out for a second-order 
optimization of UHF in Ref . [36] . Although the optimal damping factors A" and A^ may 
be formally determined by these parameters through the respective partial differentiations, 
some more consideration is necessary. 

3.2.2. Two dimensional Newton problem 

Eq. (50) can be rewritten as a second-order expansion with respect to A" and A^ 



1 



Fiti (A) = i?r' (0) + r A + -A*m, 



(56) 



where the hat is to indicate the two dimensional vectors and matrix 

A = 



A° 

A'^ 



9 



aA^A V ^ 



t 

2c^ 



(57) 
(58) 
(59) 
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The minimization of Eq. (56) leads to a simple Newton problem 

H\ = -g (60) 

as long as the Hessian is positive-definite. If this assumption is valid, the descendent A 
vector is obtained as 



The final expressions for A" and \^ thus become 

under the condition of A", A'^ e [0, 1]. Note that the neglection of t yields the essentially 
same expression as in the case of alternate ODA/UHF for each spin; this could actually 
result in a slow convergence by our experiences. When A", A'^ ^ [0, 1], the values should 
be set as unity to disable ODA. 

3.2.3. Modified Newton problem 

In the two dimensional problem of optimization, the Hessian of Eq. (59) is not re- 
stricted to be positive definite unfortunately. Namely, the two eigenvalues 



(7± = {Ca + C^) ± V(Ca - C^)2 + t^, (64) 

can take three cases (i) positive definite (0 < cr_ < cr+), (ii) saddle point (cr_ < < cr+), 
and (iii) non-positive definite ((T_ < (T+ < 0). For cases (ii) and (iii), the technique of 
shifted Hessian [29,37,38] is usable as (compare with Eq. (60)) 

{H + l^iyX^-g, (65) 
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where /i for unit matrix 1 is the shift parameter set latter. The modified solution of A are 
then obtained as 

- i2_(2c° + /x)(2c^ + Af)' ^ ^ 

Although the direction of A should be adjusted by this modification, two issues still 
remain. First, A is not the solution of Eq. (60) manifestly. Second, the length of A may 
still override the correct region of [0, 1]. These difficulties can be avoided by introducing 
a scahng relation [37] 

A = CA, Ae[0,l], (68) 
and the minimization problem of Eq. (56) is rewritten as 

<T(A(C)) = ^r'(O) + CrA + lCk*Hk C e [0, 1]. (69) 
The C is then obtained as 

C = l ^' ifAM<-rA, 

I -9*k/k*HX, otherwise. ^ ' 

The scaling factor ( can be regarded as a second damping factor consequently. Anyhow, 
the relation of A e [0, 1] as interpolation factors should be maintained. Finally, the shift 
parameter n is set as 



for < o-_ < (T+, 

(CT+ - CT_)/2 for CT_ < < CT+, (71) 
— (T_ for (7_ < (T+ < 0, 



in the actual processing. 

3.2.4. Algorithmic flow 

The concurrent ODA/UHF calculations can be described as follows. 
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Initialization: Choose an initial guess Dq and D^. Assemble Fq — F{Dq,D^), 
F§ = F{D^,D^), = J{D^) and = J{D^). Compute Eq = E{D^,D^). Set 
= = D^, ~fI = Fo«, >o = fL Jo = Jo, Jo = Jo and k ^ 0. 

Iteration: 

(a) Diagonalize F^, and assemble -Dfe+i, -D^+i via the aufbau principle. 

(b) Assemble Fock matrices F^+i and F^+i as well as Coulomb integrals Jj^_^_^ and 



J 



fc+i 



jTiQ; 



k+1 



JiDt+i), 



(c) Compute the UHF energy 

1 



Ek+i = -Tr 



?0 



(d) If the differences D^+i - and D^^^ - are enough small then go to 
termination; the energy convergence should be checked as well. 

(e) Set AZ},"+, = - D^^, and ADf+, = - D^^,. 

(f) Compute 

s'* = Tr 



Tr 



Fu^D'^u+i 

■x/3 . ■ 

Fu^Dl^i 



e = -Tr 
2 








1 

2 ' 


= ^Tr 
2 








1 

-2'- 



i = Tr 



(Jt - Jl^^ AD 



/3 

fe+1 
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(g) Compute eigenvalues: 



(h) Set a shift parameter as follows: 



for < (J- < (J+, 

(cr+ - cr_)/2 for o-_ < < 0-+, 
— cr_ for cr_ < cr+ < 0. 



(i) Set a set of tentative damping factors (A°^, A^) = (1, 1) if A" ^ [0, 1] or A'^ ^ [0, 1] 
and 



A° 



A 



{2c^ + //)s" - ts'^ 
^2- (2c° + //)(2c/3 + /x)' 

~ i2-(2c° + //)(2c/3 + //)' 



otherwise. 



(j) Set an optimal scaling factor C = 1 if A ^A < — ^*A and C = — ^*A/A ^A 
otherwise, where 



5 



and H 



2c° i 



(k) Set the fimal set of dampling factors A" = CA" and A'^ = CA^, and compute 
interpolations: 



D 
D 
F 
F 



k+l 

k+l 
a 

k+l 
k+l 



ja 



7/3 



(1 - X-)D- + A"Z^,\„ 
(1-A^)^f + A^<„ 
(1 - X'^iPl + X^F^^, + (A^ - A") [J^, - 
(1 - \^)~f[ + A^Ff^, + (A" - A'^) [j; 
(l-A")4" + AV,\i, 
(l-AOjf + A^J,V 



— 7" 

k+l "^k 
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(1) Set A; = A; + 1 and go to 2(a) 

3. Termination: Set = C^^^, = Cf+i, = D-^^, = D^^^, = F-_^„ 
= and E = Ek+i- 

As just seen, the concurrent algorithm is more comphcated than the alternate one. 

4. Test calculations 

To test the proposed ODA algorithms of both alternate and concurrent UHF cal- 
culations, we implemented them into a local version of ABINIT-MPX [39], our original 
program for the fragment molecular orbital (FMO) calculations [40] , to which UHF energy 
and its nuclear gradient had been implemented [41] as an independent work from Ref. [42]. 
We here performed the regular UHF (without FMO) calculations for four small radicals of 
the spin doublet (or single open-shell). The 6-31G* basis set [43] was used for CN, NO2 and 
(H20)3-|-OH. A hexa-aqua divalent copper complex, Cu+^-|-(H20)6, was calculated with 
the 6-31 G basis set [43], where the symmetry was imposed for the Jahn- Teller defor- 
mation due to 3d^ occupation. The geometries of four molecular systems were optimized 
by the GAUSSIAN03 program [44] at the UHF level. The extended Hiickel method was 
used to guess the initial values of AO-MO coefficients or density matrices. The convergence 
conditions of SCF iterations (cycle hmit 1000) were tightly set as \\Ek — Ek-i\\ < 10~^, 
ll^fe - Dk-i\\ < 10-s (for occupied MOs) and Max{DkSFk-i - Fk-iSD^) < 10"^ For 
testing purpose, we enforced the ODA/UHF procedure throughout (although C2-DHS 
[15] was available as well). The reference UHF energies and spin expectation values com- 
puted by GAUSSIAN03 were reproduced by ABINIT-MPX within reasonable numerical 
tolerance when converged. For comparison, the simple SCF procedure (in the sense of 
Roothaan) was adopted as well by the fixed setting A" = = 1. 

Figure 1 plots the convergence behaviors of the CN calculations. The alternate UHF 
procedure both with and without ODA shows smooth convergence, where no acceleration 
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is obtained due to a continuous resetting of the damping parameter as unity during the 
iteration. The concurrent UHF calculation without ODA fails in convergence, as expected 
from a demanding nature of this radical. In contrast, the concurrent ODA/UHF provides 
a convergence comparable to the alternate treatment. For NO2 presented in Figure 2, the 
concurrent calculation without ODA has a slow convergence, and the ODA acceleration 
works well. The behavior of the alternate calculations is similar to the case of CN. 

As seen in Figure 3, the energy lowering with the concurrent ODA procedure is rapid 
in early steps for (H20)3+OH. However, its acceleration drops off near the convergence 
in six decimal places unfortunately. This suggests that the acceleration procedure is to 
be switched to other methods such as DIIS [12,13,15,16] once certain criteria of initial 
convergence are passed. Note that the resetting of damping parameter took place for the 
concurrent calculation after the early stage. 

Figure 4 shows that both ODAs converged for Cu+^+(H20)6 but the SCF iterations 
without ODA lead to the oscillation. Notably, the concurrent version is much better, 
especially in the early stage. As denoted in the above paragraph, ODA should be switched 
to DIIS for the accelerated final convergence in production runs. We tried another initial 
guess with the diagonalization of h (core Hamiltonian) . As a result of this attempt, the 
concurrent ODA calculation converged as in the case of Hiickel guess, while the alternate 
one oscillated (data not shown). A notable merit of ODA/RHF could be a robustness 
against poor initial guesses [23], and this might be valid also for the concurrent ODA/UHF 
procedure in which the a- (3 coupling is carefully taken into account. Finally, it is a 
favorable fact that the concurrent version works better than does the alternate one for 
the four examples employed here, since the the former can be faster in processing of 
AO-integrals [25]. 

5. Summary 
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In this paper, we proposed two ODAs for UHF calculations of open-shell molecular 
systems, as extensions of the original ODA/RHF developed by Cances and Bris [23]. The 
equations associated with the Fock and density matrices were systematically derived for 
both alternate and concurrent SCF procedures. In the latter procedure, an additional 
two-dimensional Newton method was employed to determine the optimal set of damping 
factors. Test calculations were performed for four doublet radical systems. It was shown 
that the concurrent ODA has better overall performance in convergence than does the 
alternate one. This fact should be favorable since the number of integral evaluations could 
be halved in integral-direct SCF computations [25,26]. Works to fully implement the 
proposed recipes are underway for the improved performance of FMO-UHF calculations 
in the ABINIT-MPX program [39,41]. 
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Figure captions 

Figure 1. Convergence behavior of CN molecule (^E"*" state). E^. is the finally converged 
UHF energy, and Ek means the snapshot energies during the SCF iteration. The result of 
concurrent UHF calculation with ODA is labeled as "Conc/ODA" (red solid line), while 
the case without ODA is plotted as "Conc/SCF" (purple broken line). The behaviors of 
alternate UHF calculations with and without ODA are shown with labels of " Alter/ODA" 
(blue broken line) and "Alter/SCF" (green dotted line), respectively. 

Figure 2. Convergence behavior of NO2 molecule (^Ai state). The captions are the same 
as those of Figure 1. 

Figure 3. Convergence behavior of (H20)3+OH cluster (^A state). The captions are the 
same as those of Figure 1. 

Figure 4. Convergence behavior of Cu"'"^+(H20)6 complex (^Ag state). The captions are 
the same as those of Figure 1. 
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